# ------------------------------------------------------------------------------#
# Description: Create bootstrap nonlinear and number of adoption datasets
# Project:     Peer Effects in Voluntary Environmental Policies:
#              An Application to Urban Water Quality
# Author:      Daniel A. Brent | dab320@psu.edu
# Created:     2025-06-11
# ------------------------------------------------------------------------------#

# Clear environment ------------------------------------------------------------
rm(list = ls())
gc()

# Load required packages -------------------------------------------------------
library(pacman)
p_load(data.table, fixest, modelsummary, ggplot2)
options(scipen = 999)



# Load adoption data -----------------------------------------------------------
load("data/adoptions/adoptions.RData")

# Merge peer adoption buckets --------------------------------------------------
load("data/bucket/peer_adopt_buckets_any.RData")
dat <- merge(adopt, peer.adopt.bucket.any, by = c("pin", "year"))
rm(adopt, peer.adopt.bucket.any)

load("data/bucket/peer_adopt_buckets_rg.RData")
dat <- merge(dat, peer.adopt.bucket.rg, by = c("pin", "year"))
rm(peer.adopt.bucket.rg)

load("data/bucket/peer_adopt_buckets_any_rg.RData")
dat <- merge(dat, peer.adopt.bucket.any.rg, by = c("pin", "year"))
rm(peer.adopt.bucket.any.rg)

load("data/bucket/peer_adopt_buckets_cs.RData")
dat <- merge(dat, peer.adopt.bucket.cs, by = c("pin", "year"))
rm(peer.adopt.bucket.cs)

load("data/bucket/peer_adopt_buckets_both.RData")
dat <- merge(dat, peer.adopt.bucket.both, by = c("pin", "year"))
rm(peer.adopt.bucket.both)

# Merge peer eligibility buckets -----------------------------------------------
load("data/bucket/peer_elig_buckets_cs.RData")
dat <- merge(dat, peer.buckets.cs, by = c("pin", "year"))
rm(peer.buckets.cs)

load("data/bucket/peer_elig_buckets_rg.RData")
dat <- merge(dat, peer.buckets.rg, by = c("pin", "year"))
rm(peer.buckets.rg)

# Merge total peers ------------------------------------------------------------
load("data/bucket/total_peers.RData")
dat <- merge(dat, total.peers, by = "pin")
rm(total.peers)

# Merge parcel data ------------------------------------------------------------
load("data/parcel_build/parcel_build.RData")
dat <- merge(dat, parcel.build[, .(pin, zip5, sub.area, sqft, lot, beds, baths,
                                   yearbuilt, val.land, val.improve,
                                   ren.pre90, ren.90.99, ren.00.09, ren.10,
                                   water.prob)], by = "pin")
rm(parcel.build)

# Merge census data ------------------------------------------------------------
load("data/census/all_parcels_census_2010.RData")
parcels.2010[, blk := as.numeric(paste0(block, blkgrp))]
dat <- merge(dat, parcels.2010[, .(pin, tract, blkgrp, block, blk)])
rm(parcels.2010)

# Scale peer adoption variables ------------------------------------------------
dat[, peer.adopt.any.100.scale     := peer.adopt.any.100 / 100]
dat[, peer.adopt.any.rg.100.scale  := peer.adopt.any.rg.100 / 100]
dat[, peer.adopt.rg.100.scale      := peer.adopt.rg.100 / 100]
dat[, peer.adopt.cs.100.scale      := peer.adopt.cs.100 / 100]
dat[, peer.adopt.both.100.scale    := peer.adopt.both.100 / 100]

# Keep post-2010 only ----------------------------------------------------------
dat <- dat[year > 2010]

# Create average eligibility variable ------------------------------------------
dat[, elig.peers.avg.100 := elig.peers.cs.100 + elig.peers.rg.100]

# Merge shoreline data ---------------------------------------------------------
load("data/shoreline/shoreline.RData")
setnames(near.shore, "near_dist", "near_shore")
dat <- merge(dat, near.shore[, .(pin, near_shore)], by = "pin", all.x = TRUE)
rm(near.shore)

# Merge arterial road data -----------------------------------------------------
load("data/arterial/arterial.RData")
setnames(near.arterial, "near_dist", "near_road")
dat <- merge(dat, near.arterial[, .(pin, near_road)], by = "pin", all.x = TRUE)
rm(near.arterial)

# Merge mandatory GSI data -----------------------------------------------------
load("data/mandatory_gsi/peer_adopt_buckets_any_mandgsi.RData")
dat <- merge(dat, peer.adopt.bucket.any.mandgsi, by = c("pin", "year"))
rm(peer.adopt.bucket.any.mandgsi)

# Create IV Table 3 and 4 Dataset  --------------------------------------------------
table3.iv.mult <- ivDiag(data=dat[elig.cs==1 & post.rainwise==0],
                         Y = 'adopt.any',
                         D =  'peer.adopt.any.100.scale', 
                         Z = c('elig.peers.cs.100', 'elig.peers.rg.100'),
                         controls = 'peers.100',
                         FE = c('year' , 'blkgrp'),
                         cl='blkgrp',
                         bootstrap = F, run.AR = T)

table3.iv <- ivDiag(data=dat[elig.cs==1 & post.rainwise==0],
                    Y = 'adopt.any',
                    D =  'peer.adopt.any.100.scale', 
                    Z = c('elig.peers.avg.100'),
                    controls = 'peers.100',
                    FE = c('year' , 'blkgrp'),
                    cl='blkgrp',
                    bootstrap = F, run.AR = T)

table4.any.rg <- ivDiag(data=dat[elig.cs==1 & post.rainwise==0],
                        Y = 'adopt.any.rg',
                        D =  'peer.adopt.any.100.scale', 
                        Z = c('elig.peers.avg.100'),
                        controls = 'peers.100',
                        FE = c('year' , 'blkgrp'),
                        cl='blkgrp',
                        bootstrap = F, run.AR = T)

table4.cs <- ivDiag(data=dat[elig.cs==1 & post.rainwise==0],
                    Y = 'adopt.cs',
                    D =  'peer.adopt.any.100.scale', 
                    Z = c('elig.peers.avg.100'),
                    controls = 'peers.100',
                    FE = c('year' , 'blkgrp'),
                    cl='blkgrp',
                    bootstrap = F, run.AR = T)

table4.rg <- ivDiag(data=dat[elig.cs==1 & post.rainwise==0],
                    Y = 'adopt.rg',
                    D =  'peer.adopt.any.100.scale', 
                    Z = c('elig.peers.avg.100'),
                    controls = 'peers.100',
                    FE = c('year' , 'blkgrp'),
                    cl='blkgrp',
                    bootstrap = F, run.AR = T)

table4.both <- ivDiag(data=dat[elig.cs==1 & post.rainwise==0],
                      Y = 'adopt.any.rg',
                      D =  'peer.adopt.any.100.scale', 
                      Z = c('elig.peers.avg.100'),
                      controls = 'peers.100',
                      FE = c('year' , 'blkgrp'),
                      cl='blkgrp',
                      bootstrap = F, run.AR = T)

iv_diag_table3_4 <- list(table3.iv.mult,table3.iv,table4.any.rg,table4.cs,table4.rg,table4.both)
names(iv_diag_table3_4) <- c('table3.iv.mult','table3.iv','table4.any.rg',
                             'table4.cs','table4.rg','table4.both')

save(iv_diag_table3_4, file = 'output/iv/iv_diag_table3_4.RData')

rm(table3.iv.mult,table3.iv,table4.any.rg,table4.cs,table4.rg,table4.both)


# Create IV Diag Num of Peers Dataset  --------------------------------------------------

for (i in c(100,80,60,40,20)) {
  for (g in c('any','any.rg','cs','both')) {
    dat[,paste0('peer.adopt.any.',i,'.scale') := get(paste0('peer.adopt.any.',i))/100]
    dat[,paste0('elig.peers.avg',i) := (get(paste0('elig.peers.cs.',i)) + 
                                          get(paste0('elig.peers.rg.',i)))/2]
    
    assign(paste0('m.',g,'.',i), 
           ivDiag(data=dat[elig.cs==1 & post.rainwise==0],
                  Y = paste0('adopt.',g),
                  D = paste0('peer.adopt.any.',i,'.scale'),
                  Z = paste0('elig.peers.avg',i),
                  controls = paste0('peers.',i),
                  FE = c('year' , 'blkgrp'),
                  cl='blkgrp',
                  bootstrap = F, run.AR = T)
    )
  }
}

iv_diag_numpeers <- list(m.any.100,m.any.80,m.any.60,m.any.40,m.any.20,
                         m.any.rg.100,m.any.rg.80,m.any.rg.60,m.any.rg.40,m.any.rg.20,
                         m.cs.100,m.cs.80,m.cs.60,m.cs.40,m.cs.20,
                         m.both.100,m.both.80,m.both.60,m.both.40,m.both.20)
names(iv_diag_numpeers) <- c('m.any.100','m.any.80','m.any.60','m.any.40','m.any.20',
                             'm.any.rg.100','m.any.rg.80','m.any.rg.60','m.any.rg.40','m.any.rg.20',
                             'm.cs.100','m.cs.80','m.cs.60','m.cs.40','m.cs.20',
                             'm.both.100','m.both.80','m.both.60','m.both.40','m.both.20')

save(iv_diag_numpeers, file = 'output/iv/iv_diag_numpeers.RData')
rm(m.any.100,m.any.80,m.any.60,m.any.40,m.any.20,
   m.any.rg.100,m.any.rg.80,m.any.rg.60,m.any.rg.40,m.any.rg.20,
   m.cs.100,m.cs.80,m.cs.60,m.cs.40,m.cs.20,
   m.both.100,m.both.80,m.both.60,m.both.40,m.both.20)

# Create IV Diag Timing Dataset  --------------------------------------------------


# Peer Adoptions by Time
load('data/time/peer_adopt_t_any.RData')
dat <- merge(dat,peer.adopt.t.any,by=c('pin','year'))
rm(peer.adopt.t.any)

# Merge in peer adoptions (any rg)
load('data/time/peer_adopt_t_any_rg.RData')
dat <- merge(dat,peer.adopt.t.any.rg,by=c('pin','year'))
rm(peer.adopt.t.any.rg)

# Merge in peer adoptions (cs)
load('data/time/peer_adopt_t_cs.RData')
dat <- merge(dat,peer.adopt.t.cs,by=c('pin','year'))
rm(peer.adopt.t.cs)

# Merge in peer adoptions (both)
load('data/time/peer_adopt_t_both.RData')
dat <- merge(dat,peer.adopt.t.both,by=c('pin','year'))
rm(peer.adopt.t.both)

# Create IV Diag Timing Dataset  --------------------------------------------------

for (i in 1:5) {
  for (g in c('any','any.rg','cs','both')) {
    dat[,paste0('peer.adopt.any.',i,'.scale') := get(paste0('peer.adopt.any.',i))/100]
    
    assign(paste0('m.',g,'.',i),
           ivDiag(data=dat[elig.cs==1 & post.rainwise==0],
                  Y = paste0('adopt.',g),
                  D = paste0('peer.adopt.any.',i,'.scale'), 
                  Z = 'elig.peers.avg.100',
                  controls = 'peers.100',
                  FE = c('year' , 'blkgrp'),
                  cl='blkgrp',
                  bootstrap = F, run.AR = T)
    )
  }
}

iv_diag_timing <- list(m.any.1,m.any.2,m.any.3,m.any.4,m.any.5,
                       m.any.rg.1,m.any.rg.2,m.any.rg.3,m.any.rg.4,m.any.rg.5,
                       m.cs.1,m.cs.2,m.cs.3,m.cs.4,m.cs.5,
                       m.both.1,m.both.2,m.both.3,m.both.4,m.both.5)
names(iv_diag_timing) <- c('m.any.1','m.any.2','m.any.3','m.any.4','m.any.5',
                           'm.any.rg.1','m.any.rg.2','m.any.rg.3','m.any.rg.4','m.any.rg.5',
                           'm.cs.1','m.cs.2','m.cs.3','m.cs.4','m.cs.5',
                           'm.both.1','m.both.2','m.both.3','m.both.4','m.both.5')

save(iv_diag_timing, file = 'output/iv/iv_diag_timing.RData')

rm(m.any.1,m.any.2,m.any.3,m.any.4,m.any.5,
   m.any.rg.1,m.any.rg.2,m.any.rg.3,m.any.rg.4,m.any.rg.5,
   m.cs.1,m.cs.2,m.cs.3,m.cs.4,m.cs.5,
   m.both.1,m.both.2,m.both.3,m.both.4,m.both.5)


dat[, c("peer.adopt.any.5" , "peer.adopt.any.4","peer.adopt.any.3","peer.adopt.any.2",
        "peer.adopt.any.1","peer.adopt.any.rg.5","peer.adopt.any.rg.4","peer.adopt.any.rg.3",
        "peer.adopt.any.rg.2","peer.adopt.any.rg.1","peer.adopt.cs.5","peer.adopt.cs.4",
        "peer.adopt.cs.3","peer.adopt.cs.2","peer.adopt.cs.1","peer.adopt.both.5",
        "peer.adopt.both.4","peer.adopt.both.3","peer.adopt.both.2","peer.adopt.both.1",
        "peer.adopt.any.1.scale","peer.adopt.any.2.scale","peer.adopt.any.3.scale",
        "peer.adopt.any.4.scale","peer.adopt.any.5.scale") := NULL]

# Create IV Diag Distance Dataset  --------------------------------------------------

# Peer Adoptions by Distance
load('data/distance/peer_adopt_distance_any.RData')
dat <- merge(dat,peer.adopt.dist.any,by=c('pin','year'))
rm(peer.adopt.dist.any)

# Merge in peer eligibility (cs)
load('data/distance/peer_elig_distance_cs.RData')
dat <- merge(dat,peer.dist.cs,by=c('pin','year'))
rm(peer.dist.cs)

# Merge in peer eligibility (rg)
load('data/distance/peer_elig_distance_rg.RData')
dat <- merge(dat,peer.dist.rg,by=c('pin','year'))
rm(peer.dist.rg)


# Table for different radii (0.1 to 0.5 miles) 
for (i in 1:5) {
  dat[,paste0('peer.adopt.any.',i,'.scale') := get(paste0('peer.adopt.any.',i))/100]
  dat[,paste0('elig.peers.avg.',i) := (get(paste0('elig.peers.cs.',i)) + 
                                         get(paste0('elig.peers.rg.',i)))/2]
  
  assign(paste0('m',i), ivDiag(data=dat[elig.cs==1 & post.rainwise==0],
                               Y = 'adopt.any',
                               D = paste0('peer.adopt.any.',i,'.scale'), 
                               Z = paste0('elig.peers.avg.',i),
                               controls = paste0('peers.',i),
                               FE = c('year' , 'blkgrp'),
                               cl='blkgrp',
                               bootstrap = F, run.AR = T)
  )
}

iv_diag_distance <- list(m1,m2,m3,m4,m5)

names(iv_diag_distance) <- c('m1','m2','m3','m4','m5')

save(iv_diag_distance, file = 'output/iv/iv_diag_distance.RData')
rm(m1,m2,m3,m4,m5)

dat[,c("peer.adopt.any.5","peer.adopt.any.4","peer.adopt.any.3","peer.adopt.any.2",
       "peer.adopt.any.1","peer.adopt.any.1.scale","elig.peers.cs.5","elig.peers.cs.4",
       "elig.peers.cs.3","elig.peers.cs.2","elig.peers.cs.1","elig.peers.rg.5",
       "elig.peers.rg.4","elig.peers.rg.3","elig.peers.rg.2","elig.peers.rg.1",
       "peer.adopt.any.2.scale","peer.adopt.any.3.scale","peer.adopt.any.4.scale",
       "peer.adopt.any.5.scale","elig.peers.avg.1","elig.peers.avg.2",
       "elig.peers.avg.3","elig.peers.avg.4","elig.peers.avg.5") := NULL]


# Create IV Diag Table 3 & 4 (RG) Dataset  --------------------------------------------------

table3.iv.mult <- ivDiag(data=dat[elig.rg==1 & post.rainwise==0],
                         Y = 'adopt.any',
                         D =  'peer.adopt.any.100.scale', 
                         Z = c('elig.peers.cs.100', 'elig.peers.rg.100'),
                         controls = 'peers.100',
                         FE = c('year' , 'blkgrp'),
                         cl='blkgrp',
                         bootstrap = F, run.AR = T)

table3.iv <- ivDiag(data=dat[elig.rg==1 & post.rainwise==0],
                    Y = 'adopt.any',
                    D =  'peer.adopt.any.100.scale', 
                    Z = c('elig.peers.avg.100'),
                    controls = 'peers.100',
                    FE = c('year' , 'blkgrp'),
                    cl='blkgrp',
                    bootstrap = F, run.AR = T)

table4.any.rg <- ivDiag(data=dat[elig.rg==1 & post.rainwise==0],
                        Y = 'adopt.any.rg',
                        D =  'peer.adopt.any.100.scale', 
                        Z = c('elig.peers.avg.100'),
                        controls = 'peers.100',
                        FE = c('year' , 'blkgrp'),
                        cl='blkgrp',
                        bootstrap = F, run.AR = T)

table4.cs <- ivDiag(data=dat[elig.rg==1 & post.rainwise==0],
                    Y = 'adopt.cs',
                    D =  'peer.adopt.any.100.scale', 
                    Z = c('elig.peers.avg.100'),
                    controls = 'peers.100',
                    FE = c('year' , 'blkgrp'),
                    cl='blkgrp',
                    bootstrap = F, run.AR = T)

table4.rg <- ivDiag(data=dat[elig.rg==1 & post.rainwise==0],
                    Y = 'adopt.rg',
                    D =  'peer.adopt.any.100.scale', 
                    Z = c('elig.peers.avg.100'),
                    controls = 'peers.100',
                    FE = c('year' , 'blkgrp'),
                    cl='blkgrp',
                    bootstrap = F, run.AR = T)

table4.both <- ivDiag(data=dat[elig.rg==1 & post.rainwise==0],
                      Y = 'adopt.any.rg',
                      D =  'peer.adopt.any.100.scale', 
                      Z = c('elig.peers.avg.100'),
                      controls = 'peers.100',
                      FE = c('year' , 'blkgrp'),
                      cl='blkgrp',
                      bootstrap = F, run.AR = T)

iv_diag_table3_4_rg <- list(table3.iv.mult,table3.iv,table4.any.rg,table4.cs,table4.rg,table4.both)
names(iv_diag_table3_4_rg) <- c('table3.iv.mult','table3.iv','table4.any.rg',
                                'table4.cs','table4.rg','table4.both')
save(iv_diag_table3_4_rg, file = 'output/iv/iv_diag_table3_4_rg.RData')
rm(table3.iv.mult,table3.iv,table4.any.rg,table4.cs,table4.rg,table4.both)

# Create IV Diag Basins Dataset  --------------------------------------------------

# Base
m1 <- ivDiag(data=dat[elig.cs==1 & post.rainwise==0],
             Y = 'adopt.any',
             D =  'peer.adopt.any.100.scale', 
             Z = c('elig.peers.avg.100'),
             controls = c('peers.100'),
             FE = c('year' , 'blkgrp'),
             cl='blkgrp',
             bootstrap = F, run.AR = T)


# Parcel controls
m2 <- ivDiag(data=dat[elig.cs==1 & post.rainwise==0],
             Y = 'adopt.any',
             D =  'peer.adopt.any.100.scale', 
             Z = c('elig.peers.avg.100'),
             controls = c('peers.100','sqft', 'lot','beds','baths','yearbuilt',
                          'ren.pre90','ren.90.99','ren.00.09','ren.10', 'water.prob'),
             FE = c('year' , 'blkgrp'),
             cl='blkgrp',
             bootstrap = F, run.AR = T)

# Shore control
m3 <- ivDiag(data=dat[elig.cs==1 & post.rainwise==0],
             Y = 'adopt.any',
             D =  'peer.adopt.any.100.scale', 
             Z = c('elig.peers.avg.100'),
             controls = c('peers.100','sqft', 'lot','beds','baths','yearbuilt',
                          'ren.pre90','ren.90.99','ren.00.09','ren.10', 
                          'water.prob','near_shore'),
             FE = c('year' , 'blkgrp'),
             cl='blkgrp',
             bootstrap = F, run.AR = T)


# Road control
dat[,close.road := ifelse(!is.na(near_road),1,0)]
dat[is.na(near_road),near_road := 0]
dat[,close.road.near.road := close.road*near_road]
m4 <- ivDiag(data=dat[elig.cs==1 & post.rainwise==0],
             Y = 'adopt.any',
             D =  'peer.adopt.any.100.scale', 
             Z = c('elig.peers.avg.100'),
             controls = c('peers.100','sqft', 'lot','beds','baths','yearbuilt',
                          'ren.pre90','ren.90.99','ren.00.09','ren.10', 
                          'water.prob','close.road','close.road.near.road'),
             FE = c('year' , 'blkgrp'),
             cl='blkgrp',
             bootstrap = F, run.AR = T)

# Mandatory GSI
m5 <- ivDiag(data=dat[elig.cs==1 & post.rainwise==0],
             Y = 'adopt.any',
             D =  'peer.adopt.any.100.scale', 
             Z = c('elig.peers.avg.100'),
             controls = c('peers.100','sqft', 'lot','beds','baths','yearbuilt',
                          'ren.pre90','ren.90.99','ren.00.09','ren.10', 
                          'water.prob','peer.adopt.any.mandgsi.100'),
             FE = c('year' , 'blkgrp'),
             cl='blkgrp',
             bootstrap = F, run.AR = T)

# DEM
m6 <- ivDiag(data=dat[elig.cs==1 & post.rainwise==0],
             Y = 'adopt.any',
             D =  'peer.adopt.any.100.scale', 
             Z = c('elig.peers.avg.100'),
             controls = c('peers.100','sqft', 'lot','beds','baths','yearbuilt',
                          'ren.pre90','ren.90.99','ren.00.09','ren.10', 
                          'water.prob','elevation_ft'),
             FE = c('year' , 'blkgrp'),
             cl='blkgrp',
             bootstrap = F, run.AR = T)

# All
m7 <- ivDiag(data=dat[elig.cs==1 & post.rainwise==0],
             Y = 'adopt.any',
             D =  'peer.adopt.any.100.scale', 
             Z = c('elig.peers.avg.100'),
             controls = c('peers.100','sqft', 'lot','beds','baths','yearbuilt',
                          'ren.pre90','ren.90.99','ren.00.09','ren.10', 
                          'water.prob','near_shore' ,'near_road',
                          'peer.adopt.any.mandgsi.100','elevation_ft'),
             FE = c('year' , 'blkgrp'),
             cl='blkgrp',
             bootstrap = F, run.AR = T)

# Drop shore
dat[is.na(near_shore)]
quantile(dat$near_shore,0.2)
m8 <- ivDiag(data=dat[elig.cs==1 & post.rainwise==0 & near_shore >= quantile(dat$near_shore,0.2)],
             Y = 'adopt.any',
             D =  'peer.adopt.any.100.scale', 
             Z = c('elig.peers.avg.100'),
             controls = c('peers.100','sqft', 'lot','beds','baths','yearbuilt',
                          'ren.pre90','ren.90.99','ren.00.09','ren.10', 
                          'water.prob','near_shore' ,'near_road',
                          'peer.adopt.any.mandgsi.100','elevation_ft'),
             FE = c('year' , 'blkgrp'),
             cl='blkgrp',
             bootstrap = F, run.AR = T)
# Drop road
dat[is.na(near_road)]
quantile(dat$near_road,0.2)
m9 <- ivDiag(data=dat[elig.cs==1 & post.rainwise==0 & near_shore >= quantile(dat$near_road,0.2)],
             Y = 'adopt.any',
             D =  'peer.adopt.any.100.scale', 
             Z = c('elig.peers.avg.100'),
             controls = c('peers.100','sqft', 'lot','beds','baths','yearbuilt',
                          'ren.pre90','ren.90.99','ren.00.09','ren.10', 
                          'water.prob','near_shore' ,'near_road',
                          'peer.adopt.any.mandgsi.100','elevation_ft'),
             FE = c('year' , 'blkgrp'),
             cl='blkgrp',
             bootstrap = F, run.AR = T)

iv_diag_robust <- list(m1,m2,m3,m4,m5,m6,m7,m8,m9)
save(iv_diag_robust, file = 'output/iv/iv_diag_robust.RData')
rm(m1,m2,m3,m4,m5,m6,m7,m8,m9)

# Blkgrp
m1 <- ivDiag(data=dat[elig.cs==1 & post.rainwise==0],
             Y = 'adopt.any',
             D =  'peer.adopt.any.100.scale', 
             Z = c('elig.peers.avg.100'),
             controls = c('peers.100','sqft', 'lot','beds','baths','yearbuilt',
                          'ren.pre90','ren.90.99','ren.00.09','ren.10', 'water.prob'),
             FE = c('year' , 'blkgrp'),
             cl='blkgrp',
             bootstrap = F, run.AR = T)
# Basin
m2 <- ivDiag(data=dat[elig.cs==1 & post.rainwise==0],
             Y = 'adopt.any',
             D =  'peer.adopt.any.100.scale', 
             Z = c('elig.peers.avg.100'),
             controls = c('peers.100','sqft', 'lot','beds','baths','yearbuilt',
                          'ren.pre90','ren.90.99','ren.00.09','ren.10', 'water.prob'),
             FE = c('year' , 'basin'),
             cl='basin',
             bootstrap = F, run.AR = T)

# Zip5
m3 <- ivDiag(data=dat[elig.cs==1 & post.rainwise==0],
             Y = 'adopt.any',
             D =  'peer.adopt.any.100.scale', 
             Z = c('elig.peers.avg.100'),
             controls = c('peers.100','sqft', 'lot','beds','baths','yearbuilt',
                          'ren.pre90','ren.90.99','ren.00.09','ren.10', 'water.prob'),
             FE = c('year' , 'basin'),
             cl='basin',
             bootstrap = F, run.AR = T)

# Sub area
m4 <- ivDiag(data=dat[elig.cs==1 & post.rainwise==0],
             Y = 'adopt.any',
             D =  'peer.adopt.any.100.scale', 
             Z = c('elig.peers.avg.100'),
             controls = c('peers.100','sqft', 'lot','beds','baths','yearbuilt',
                          'ren.pre90','ren.90.99','ren.00.09','ren.10', 'water.prob'),
             FE = c('year' , 'sub.area'),
             cl='sub.area',
             bootstrap = F, run.AR = T)

# Tract
m5 <- ivDiag(data=dat[elig.cs==1 & post.rainwise==0],
             Y = 'adopt.any',
             D =  'peer.adopt.any.100.scale', 
             Z = c('elig.peers.avg.100'),
             controls = c('peers.100','sqft', 'lot','beds','baths','yearbuilt',
                          'ren.pre90','ren.90.99','ren.00.09','ren.10', 'water.prob'),
             FE = c('year' , 'tract'),
             cl='tract',
             bootstrap = F, run.AR = T)


# Basin-year
class(dat$blkgrp)
class(dat$year)
dat[,basin.year := paste0(as.character(basin),as.character(year))]
m6 <- ivDiag(data=dat[elig.cs==1 & post.rainwise==0],
             Y = 'adopt.any',
             D =  'peer.adopt.any.100.scale', 
             Z = c('elig.peers.avg.100'),
             controls = c('peers.100','sqft', 'lot','beds','baths','yearbuilt',
                          'ren.pre90','ren.90.99','ren.00.09','ren.10', 'water.prob'),
             FE = c('basin.year'),
             cl='basin',
             bootstrap = F, run.AR = T)

# Blkgrp-year
class(dat$blkgrp)
class(dat$year)
dat[,blkgrp.year := paste0(as.character(blkgrp),as.character(year))]
m7 <- ivDiag(data=dat[elig.cs==1 & post.rainwise==0],
             Y = 'adopt.any',
             D =  'peer.adopt.any.100.scale', 
             Z = c('elig.peers.avg.100'),
             controls = c('peers.100','sqft', 'lot','beds','baths','yearbuilt',
                          'ren.pre90','ren.90.99','ren.00.09','ren.10', 'water.prob'),
             FE = c('blkgrp.year'),
             cl='blkgrp',
             bootstrap = F, run.AR = T)

iv_diag_spatial <- list(m1,m2,m3,m4,m5,m6,m7)
save(iv_diag_spatial, file = 'output/iv/iv_diag_spatial.RData')
rm(m1,m2,m3,m4,m5,m6,m7)


# Identify basins with large and negative effects
large <- c('E. Ballard/W. Phinney Ridge','Madison Park','Queen Anne','Fremont/Wallingford')
neg <- c('Portage Bay','Duwamish','Barton/Fauntleroy')
## not %in%
'%!in%' <- function(x,y)!('%in%'(x,y))  ## create your own function

m1 <- ivDiag(data=dat[elig.cs==1 & post.rainwise==0],
             Y = 'adopt.any',
             D =  'peer.adopt.any.100.scale', 
             Z = c('elig.peers.avg.100'),
             controls = c('peers.100'),
             FE = c('year' , 'blkgrp'),
             cl='blkgrp',
             bootstrap = F, run.AR = T)

m2 <- ivDiag(data=dat[elig.cs==1 & post.rainwise==0 & basin %!in% large],
             Y = 'adopt.any',
             D =  'peer.adopt.any.100.scale', 
             Z = c('elig.peers.avg.100'),
             controls = c('peers.100'),
             FE = c('year' , 'blkgrp'),
             cl='blkgrp',
             bootstrap = F, run.AR = T)

m3 <- ivDiag(data=dat[elig.cs==1 & post.rainwise==0 & basin %!in% neg],
             Y = 'adopt.any',
             D =  'peer.adopt.any.100.scale', 
             Z = c('elig.peers.avg.100'),
             controls = c('peers.100'),
             FE = c('year' , 'blkgrp'),
             cl='blkgrp',
             bootstrap = F, run.AR = T)

m4 <- ivDiag(data=dat[elig.cs==1 & post.rainwise==0 & basin %!in% c(large,neg)],
             Y = 'adopt.any',
             D =  'peer.adopt.any.100.scale', 
             Z = c('elig.peers.avg.100'),
             controls = c('peers.100'),
             FE = c('year' , 'blkgrp'),
             cl='blkgrp',
             bootstrap = F, run.AR = T)


iv_diag_basins <- list(m1,m2,m3,m4)
save(iv_diag_basins, file = 'output/iv/iv_diag_basins.RData')
